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Abstract. We propose and analyze an algorithm for the solution of the L^-subgradient flow of 
the total variation functional. The algorithm involves no regularization, thus the numerical solution 
preserves the main features that motivate practitioners to consider this type of energy. We propose 
an iterative scheme for the solution of the arising problems, show that the iterations converge, and 
develop a stopping criterion for them. We present numerical experiments which illustrate the power 
of the method, explore the solution behavior, and compare with regularized flows. 
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1. Introduction. This work is concerned with the approximation of solutions 
to the following (formal) initial boundary value problem; 



Here fl is an open, bounded and connected subset of W'' with d > 1; dil denotes the 
boundary of and n its exterior unit normal; and T > is a positive and finite 
time. This equation and related ones are commonly known as very singular diffusion 
equations (see [5S1[57]) since in flat regions, i.e., |Vu| = 0, the diffusion is so strong 
that it is not a local effect anymore. 

Beginning with the seminal paper j34) equations of this class have received con- 
siderable attention from the image processing community, since such models preserve 
discontinuities while removing noise and other artifacts (see also pTl IT51 [TBI [T7[ f5T] V 
In addition, such equations appear in the modeling of grain boundary motion |29j : 
facet formation and evolution [23 ; electroniigration |25j and various other problems 
stemming from materials science. This clearly shows that the development of efficient 
and accurate numerical schemes for the solution of this class of problems is of ex- 
treme importance. To the best of our knowledge however, the techniques advocated 
in the literature for the solution of these equations as a rule involve a regularization 
somewhat related with replacing the singular term by 



see, for instance, [9l [HI [211 [22l |30] . The disadvantages of this approach are twofold: 
first, although these methods have been shown to converge, there is no clear under- 
standing of the relation between the regularization parameter e and the discretization 
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(1.1) 



|Vu|, = ye2 + |Vu|2, e>0; 



(1.2) 
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parameters; second, the regularization destroys certain fundamental features of the 
solutions which motivate the introduction of such models in the first place. 

A notable exception to the trend mentioned above is the work [7] which, inspired 
by the ideas advanced in |12j . develops a finite element scheme for total variation 
minimization which involves no regularization. In this work we adapt and extend the 
ideas presented in [7| to the study of total variation flows. We propose and analyze 
an unconditionally stable and convergent discretization scheme for the approximation 



of (1.1). In addition, we study an iterative scheme for the solution of the discrete 
problems and develop an a posteriori error estimator that provides a robust stop- 
ping criterion for the iterative scheme that guarantees that, although we only have 
approximate solutions, the convergence properties of the method are not affected. 

This work is organized as follows. The notation and conventions are set in Sec- 
tion [2j In j ]2.1| we recall the definition and main properties of functions of bounded 
variation and, in addition, we present an approximation result for these functions. 
This will be our workhorse during the derivation of error estimates. The proper frame- 



work to understand (1.1), i.e., total variation fiow and its properties are described in 
Section [3) Time and space discretization are discussed in Section |4j where we begin 
by reviewing results on implicit semidiscretizations of gradient flows in Hilbert spaces. 
Then we provide a general theory for fully discrete subgradient flows in Hilbert spaces, 
which we later apply to our problem. One of the salient novelties of this general the- 
ory is that we allow for modifications in the energy which can be used, for instance, 
to take into account the effects of quadrature. In Section [5] we propose an iterative 
scheme for the solution of the fully discrete problems, and devise a stopping criterion 
for the iterations, which guarantees that the convergence orders are not damaged. 
Section [6] contains a series of numerical experiments, which illustrate and extend the 
properties and theory for the developed scheme. Finally, in Section [7j we apply the 
approximation result of §2.1| to the problem of total variation minimization and prove 
convergence estimates similar to those available in the literature, but with reduced 
regularity assumptions. 

2. Notation and Preliminaries. We denote by J7 a bounded domain in 
with d>l. The boundary of Vl is denoted by dO. and we assume that G C°'^. We 
set T > to be a finite final time. As usual, we denote by LP(n) the space of Lebesgue 
integrable functions with exponent p S [1, oo] and by ($7), s S K, the usual Sobolev 
spaces. Spaces of vector valued functions and their elements will be represented with 
boldface characters. Recall that the following interpolation inequality holds (cf. |24j ) 

\MLr <\\w\\l,\\w\\\-\ p,gG[l,oo], i = ^ + ^, sG[0,l]. (2.1) 

Whenever i? is a normed space, we denote its norm by || • \\e and its dual by 
E' . The L^-inner product will be denoted by (•,•). Function spaces of vector- valued 
functions will be denoted by boldface characters. For function spaces, if it is clear 
from the context, we will omit the domain of definition. We denote by *Bi(i?) the 
unit ball in E, i.e., the set ^i{E) = {x € E : \\x\\e < 1} . Let : [0,T] £; be a 
measurable function in the Bochner sense, then, for p G [l,oo], we define 

UVlhe)"^ ( li'/'WIlfidi, p < oo, |l(/)||ioc(£;) = ess sup|l(/)(i)||is. 
JO te[o,T] 

We will denote by the time derivative of 4>. 
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To deal with time discretization, we introduce a time-step At > 0, for simplicity 
assumed constant. Then we partition the time interval via t'' = kAt with k = 0,K 
and K = \T/At\. We use the notation 0^* = |(/)'^}^_^, introduce the time increment 
operator 

50*^- = (2.2) 

and the extrapolation operator 

^*,k+i ^ ^ ^ ^*,fc+i ^ j^g^k^ j^^Q (2.3) 

For 0^* C E and p G [1, cxd], we introduce the (semi)norms 

k=0 k=l 

with the usual modification for p = oo. Given a sequence 0^* C E we will want to 
be able to compare it to functions defined on [0,T]. To this end, we define the Rothe 
interpolant (f>, that is, the piecewise linear function 

?(i) = ^0'+' + (2.4) 

Notice that, by construction, (f) £ C^'^{E), and = Recall also 

the well-known summation by parts formula: for every c L^(f2), 



E {{^'^H'''') + {s<i>'^\^'+')) = (<^^, - r) ■ (2.5) 

fe=0 

We will carry out the space discretization with finite element techniques. In other 
words, given we introduce a so-called triangulation Th = {T} as a collection of cells 
that satisfy the usual conformity and shape regularity assumptions [HI [^U] , and are 
such that Cl = Utgt; parametrize our collection of triangulations via h = 

max{diam(T) : T € Th}- For simplicity, we assume that each T is the isoparametric 
image of a so-called reference cell, which can be either T — [—1, 1]'', in which case 
we call the cells cubic; or T = |(a;i, . . . , Xd) G M'* : a^i > X]f=i < l| , which are 

called simpHces. We denote by {zj,t} the vertices of the cell T and set JVt = #{zj.T}- 
Clearly, Afx = 2^^ for cubes and At = d 4- 1 for simplices. 
We define the finite element space 

Wh = {wh e c%n) : vh\T er}ci w^{n), (2.6) 

where P = Qi for cubes and V — Pi for simplices. Here Qi denotes the space of 
polynomials of degree at most one in each variable and Pi the space of polynomials 
of total degree not greater than one. 

For a function w such that w\t G C'^(T'), we define its local Lagrange interpolant 
IhW by 



Ihw\T : Ihwlrizj^r) = w\t{zj,t), j = 1,7Vt, WeTh- 



(2.7) 
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This operator satisfies 

\\w ~Ihw\\Lp(T)+h\\^(.w ~Ihw)\\Lp{T) < ch'^\\D'^w\\Lp(T), (2.8) 

see, e.g. [2], for a proof. We remark also that if w S C°(il), then XhW G Yh- It will be 
also necessary to introduce the Clement interpolant (cf. [TS]) Ilh ■ L^{fl) — > V/j. This 



operator enjoys approximation properties similar to (2.8), the only difference being 
that the domain on the right hand side is a neighborhood of T. More importantly, 
the operator is stable under any Sobolev norm, i.e., 

WUhwWw^ <c\\w\\w^, s>0, [l,oo]. (2.9) 

We will denote by c a constant whose value might change at each occurrence. 

2.1. Functions of Bounded Variation and their Approximation. We say 

that a function w E L^{il) belongs to the space BV{n) (is of bounded variation) if 
its derivative, Dw, in the sense of distributions is a Radon measure. In other words, 
|Dw|(r2) < oo, where for any Borel set A C fl, 

\Bw\{A)= Slip I^J^wVq: qeC^{A), ||q||L~(A) < 1 

The space BV{ft) endowed with the norm ||w||_By = II'^'IIli + |Dw|(r2) is a Banach 
space. For more details on this space, we refer to [TJ|3H]. Let us present a result on 
approximation by smooth functions, in the case of a star-shaped domain. 

Proposition 2.1 (Approximation of BV functions). Assume that H. is bounded, 
star-shaped with respect to a point and dfl G C^'^ . If w E BV{^), then for every e > 
there exists a C°°(il) such that 



<e\Dw\{n), \\Ww,\\li < {l + ce)\Dw\{n), \\D'^w,\\li < ce-^\Dw\{n). 



Proof Without loss of generality, we can assume that G and that fl is 
star-shaped with respect to 0. For e > we define 

{y eR"^ -.y ^ {1 + e)x, x e Q} , 

and notice that and f2e are related via a bijective and Lipschitz, in fact linear, 
transformation with Lipschitz constant 1 + e and Jacobian (1 + eY ^ 1 + ce. For 
w e BV{il) we define € BV{Qf^) via Vf_{y) — w{j^) and, for x & V,, w^{x) = 
Ve * Pe{x), where is a smooth convolution kernel such that, for ever \ < p < oo, 
WVpeWhp < ce-^i+'^/p') withp' =p/{p-l). 

If w e C^(fi) then, clearly, \\w - w^Wl^ < f]\^w\\Li = e|Dw|(fi) and 

\\Vw,\\li < \Bv,\{n,) <{l + ce)\\Vw\\Li = (1 + ce)\Bw\{n). 

We now recall that smooth functions are dense in BV{ft) under strict convergence, 
(cf. [DISH]). In other words, given w £ BV{fl) there is a sequence {wn}neN C C°°{fl) 
such that 

lim \\w — Wn\\L^ — 0, hmsup iDuinK^l) < |Du;|(r2). 

n— ^oo — 

Applying the argument given above to elements of this sequence and then passing to 
the limit we obtain the first two inequalities. 
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We use Young's inequality for convolutions [24l Proposition 8.7] to obtain 
< ||Vu;e|lLi||Vp,|Ui < c^^^\Bw\{n) < -\Bw\{n), 



which concludes the proof. □ 

Remark 2.2 (Approximation in L^-spaces). Observe that, in the setting of 
Proposition \2. 1\ Young's inequality also implies 

PVIIlp < llVwellLillVp.lliP < ce-(l+'^/p')|D^i;|(^!), 

for any p G [1, oo] with p' — p/{p — 1)- In the sequel, however, we shall avoid using 
this bound since it would lead to d-dependent error estimates. 

3. The Total Variation Flow. Let us define the functional : L'^ (n) ^ R by 

f\Dwm), weL'mnBVin), 

wilt;) = < „ (3.1) 

\+oo, w e L^{n)\BV{n), 

It is not difficult to show that is convex and lower semicontinuous. Then, one can 
define the subdiffercntial of (see [2, 6, 8, lOl UHl |26l [27]) and study its subgradient 
flow, i.e., we seek for a function u : [0,T] — > L^(ri) such that 

ut e -9*(u), (3.2) 

or, equivalently, 

(uf,u-w)+*(u)-*(w) <0, yw(zL^{n). (3.3) 



It is in this sense that (1.1) is going to be understood and analyzed. Existence of 
solutions to (3.2) can be obtained with the help of the theory of maximal monotone 
operators, see [51 [TO]. 

Remark 3.1 (Dirichlet boundary conditions). The definition that we have pro- 



vided corresponds to imposing Neumann boundary conditions as in (1.1 1. The issue 



of how to impose Dirichlet boundary conditions is a delicate one since the trace of 



a BV{n) function is in L^{dfl) (cf. J71 \3^). In addition, the functional (3.1) has 
linear growth, so that the only possible way to impose Dirichlet boundary conditions is 
with the relaxed energy ^'(w) + \w~ g\; see fB \28^ for details. The introduction of 
this additional non- differentiate term greatly complicates the analysis. However, in 
the particular situation when £7 is convex, the boundary data is time independent and 
continuous, i.e., g{x,t) — g{x) e C''\dQ.), and the initial data is compatible with the 
boundary data in the sense that Uo|ao = 9, then the solution u to the subgradient flow 
with the relaxed energy satisfies u(-,t)|(3si = 9[')j for all t € (0,T]; see \2iA Lemma 
4-.1]. This can be realized, for instance, in the case of homogeneous Dirichlet boundary 
conditions (g = 0) and compactly supported initial data. Under this particular setting 
all the results we present will also hold. 

Reference [21 Theorem 2.16] shows that problem ( |3.2[ ) possesses a L'^~-L°° regu- 
larizing effect, that is if Uq € L''(ri), then ii{t) G L°°{n) for all t > 0. Let us show 
that solutions to this problem also satisfy a maximum principle. 

Theorem 3.2 (Maximum principle for TV flow). Assume that uq e L°°{rt), then 



u, solution of is such that u G L°°{[0,T],L°°{n)) and 
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Proof. Since Uq € L°°{n) we can define Uq = ess sup^gj^ Uo(a;), and w{t) = 
u(i) VUq. Since suppDw = {x £ V, : u > Uq} and, on this set, DiD — Du, we can 
conclude that w e L^{fl)nBV{n) and \E'(w) < *(u). Setting w = win ([33| we obtain 



(ut,OV (u —Uq)) < 0, which, since is constant, imphes — Z^o) V 0||^2 < 0. 

Given that (uq — Uq) V = this imphes the resuh. □ 

4. Discretization. In this section we introduce and analyze a fully discrete 



scheme for the approximation of solutions to (3.3). We begin with a semidiscrete 
(continuous in space and discrete in time) scheme for (3.3 1 which can then be analyzed 
using standard results from the literature (cf. ^2 [HQ). Then we develop a theory for 
fully discrete subgradient flows in Hilbert spaces and discuss the effect of introducing 
a discrete energy and a perturbation on the right hand side. We will provide sufficient 
compatibility conditions between the space discretization and the discrete energy to 
guarantee convergence. These results constitute a general and, as far as we know, novel 
approach to the study of fully discrete schemes for subgradient flows and evolution 
variational inequalities. The main application of these results will be, of course, a 



fully discrete scheme for (3.3) 



4.1. A Semidiscrete Scheme for TV Flows. We introduce a sequence {m'^*} 
contained in L^(i7) n BV{il) with u'-^ = uq that solves: 

^^^,u'^+i-it;^+*(u'=+i)-*(«;)<0 yweL^in). (4.1) 

Existence and uniqueness is guaranteed by the convexity and lower semicontinuity of 
^; see [HOITnillH]. A priori estimates, as well as a maximum principle are established 
in the next result. 

Proposition 4.1 (Semidiscrete stability). Let {u^*} solve ( |4.1[ ). Ifu° e L'^{^), 
then 

K 

||u^*||,^.o(^.) + + AtY^\^u\n) < c\\uYl^. (4.2) 

fe=i 

IfvP e L'^{n) n BViVt), then the flow is monotone, i.e., 5'(m'=+1) < *(u'=) < *(m") 
for all k > 0; and, moreover, 

+ Atvl/(u^') < At^iu°). (4.3) 
Ifu° e L°°{n) n BV{n), then u^^ C L°°{n) and 



Proof. To obtain (4.2) it suffices to set w = on (4.1) and add over k. To obtain 



(4.3) we set w = u and add over k. The maximum principle is obtained mutatis 



mutandis the proof of Theorem 3.2 □ 



The convergence properties of (4.1 ) are a consequence of standard results [331135 



For convenience we summarize them below. 

Corollary 4.2 (Convergence of semidiscrete TV flow). 
L^{n)r]BV{n), then 



that 



At*(u^) < c(||un 



Uq 



2 

L2 
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where u is the piecewise linear function defined in (2.4) and a — 1. If, in addition, 
9*(uo) n L^in) ^ 0, then a = 2. 

Remar k 4.3 (Error estimates for semidiscrete flows) . Notice that the conclusion 



of Corollary 4-2 gives an error estimate of order 0{At^^^) under the sole assumption 
that the initial condition has finite energy. Under this regularity the result is optimal, 
see '35\ Theorem 5]. Observe also that the assumption d'^ (uq) Ci (fl) 7^ 0, necessary 
to obtain a 0{At) estimate, is not a regularity but rather a compatibility assumption. 
We refer the reader to 121 for a partial characterization of the subdifferential of the 
total variation. 

4.2. Fully Discrete Schemes for Subgradient Flows in Hilbert Spaces. 

Here we present and analyze of a fully discrete implicit Eulcr method for subgradient 
flows in Hilbert spaces. The main novelty of our approach is the treatment of the 
space discretization and that we allow for perturbations of the energy, as well as of 
the right hand side. 

Let _ff be a Hilbert space with inner product (•,•), which induces the norm || • \\h. 
Assume that the functional ^ : 'D{^) C if — )■ M is convex, lower semicontinuous and 
I?(S') = H. Then, for every w G the subdifferential d^{w) C H is not empty. 

For details we refer to [6l [TOl [18]. We want to study its subgradient flow: Find 
V : [0, T] — > with v|t=o = Vq, such that 

vt + d^iv) 3 0, (4.4) 

or, equivalently, 

(vt,v-u;) +S^(v) -S^(w) < 0, yvj€H. (4.5) 



The existence and uniqueness of a solution to (4.4 1 or (4.5 1 follows the theory of 
maximal monotone operators, [6l llOj. 

To discretize in time we consider the Euler method. In other words, we search 
for v^* C iJ, w° = Vq, such that 



V 



At 



''+^ -w)+^{v''+^)-d{w)<0, yweH. (4.6) 



The theory of [331 provides, under the assumption that v'^ G a 0{At 



1/2N 



error estimate as in Corollary 4.2 



We now introduce the space discretization. Let {"H/ilh^o be a family of (finite 
dimensional) subspaces of H. To be able to handle perturbations on the energy 
induced by the spatial discretization we assume that, for each h > 0, we have a 
convex and lower semicontinuous functional ■ T~Lh — >■ K that is monotone, in the 
sense that 

^h{wh) > Uwh), e Hh- (4.7) 

Moreover, we assume that the spaces possess suitable approximation properties. 
In other words, there is a dense subspace W ^ H, an operator Ch ■ W ^ Hh and 
functions Ei & C ([0,oo), [0,cx))), £^(0) = 0, z = 1,2, such that 

\\ChW^w\\H <ei{h)\\w\\w, VweW, (4.8) 

and this approximation is asymptotically energy diminishing, i.e., 

dhiChw) - d{w) < e2ih)\\w\\w, yw e W. (4.9) 
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We consider the following fully discrete problem: Find w^* C Tih such that 
/ /iv''+^ \ 



At 



(4.10) 



for all Wh € Hh, where p C -ff is a perturbation. 

We begin our analysis of the fully discrete method (4.10) with an a priori bound 
on the increments of the sequence u^*. 



Lemma 4.4 (Stability of derivatives). The solution w^* to (4.101 satisfies 



Proof. Set Wfi = v\ in (4.10). Multiplying by Ai, adding over k, and using 
Young's inequality, we obtain the result. □ 

The sequence p'^* is meant to be a perturbation induced by either discretization 
or the solution procedure. For this reason we shall assume that 



„At| 



IIP 



HH) < cAt 



1/2 



(4.11) 



Based on this estimate, the error analysis proceeds as follows. 

Theorem 4.5 (A priori error analysis). Let w"^* be the solution to (4.6) and v^^ 



the solution of ( |4.10 |. Assume z;^* £ i°°{W). If the operator Ch satisfies (4.8 1 and 
(4.9); the discrete energies satisfy ( |4.7[ ); S'/i(u°) < +oo uniformly in h; and the 
perturbations p^* satisfy (4.111, then there exists a constant c > proportional to T 
such that 

\\v^' - v^j^^H) < - vIWh + c(ei(M + S2ih))\\v^'\U^^w) + cAt. 



Proof. To simplify notation, as usual, we will denote e'^' = — v\. Set w = v^^^ 
in (4.6) and Wh — ChV^^^ in (4.10) and add the results. Using the monotonicity 
property (4.7) we obtain 



(^^,e'+'^<UC,.v'+')-Siv'^+')-{p'^+\e'^+') 



1 8v^+'^ 



\ 



Using Lemma 4.4 and (4.11) we obtain At ^ X^aLi ll'^^^lllf — whence (4.8) and 
yield 

(5e'=+i,e^-+i) <cAt[£2(/^) + £i(/i)]^'=+^lk + A<||/+i||H||e'+i|k. 

Since the sequence ||e*^*||^f is finite we can assume that its maximum is attained for 
k — K. Summing the above inequality over A: = 1, k — 1 we deduce 

K 

fe=i 



An application of Young's inequality, together with assumption (4.11) allows us to 
conclude. □ 
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4.3. Fully Discrete Scheme for TV Flows. Let us specialize the ideas pre- 
sented in |4.2|to the case of TV flows. To do so, we will work on the discrete spaces 



(2.6), define discrete energies and approximation operators, and verify that assump- 



tions (4.7)-(4.9l are satisfied. The conclusion of Theorem 4.5 will then allow us to 



obtain error estimates. 

In this setting, we consider the subgradient flow: Find u^* C that solves 



At 



(4.12) 



The discrete energy is defined by 



*h(w/i) = 



Ten 



(4.13) 



where the operator Xh was defined in (2.7 1. Notice that we have effectively replaced 



the total variation seminorm by a quadrature formula. Indeed, we can rewrite as 



'T^'T A 1 



WjlVu)ftl(zj-T), 



T<en j=i 



where Aj.t is the coordinate basis function associated with node z^j- and we denote 
the weights by Wj = \r\~^ \j,\jj:'{z)Az. Notice that this modification fits into the 



framework described in |4.2| Its utility will become clear in the following paragraph. 
The approximation operator is defined as 



(4.14) 



where 11^ is the Clem ent interpolation operator and denotes a regularization of w 
as in Proposition 2.1 with e = Y?l^ . Proving error estimates reduces to verification of 



with e 
the hypotheses of §4!2| 

Corollary 4.6 (Convergence of TV flow). Let u^* be the solution of (|4.1|) and 



the solution of (4.12) with the discrete energy given in (4.13). If Q is star shaped 



with respect to a point; u'^ € L°°{fl) H BV{rt); and ^'^(mJ'j) < c < +oo uniformly in 
h, then 

\\u^'^utre^^,.)<\\n'-v^Mh+ch'/'. 



Proof. Notice, first of all, that the given assumptions on u° translate into u'^* € 
£'^{L°°{n) n BVifl)). Consequently, we set W = BV[n) n L°^{n) and 

\\w\\w = IIw'IIl^" + |Dw|(il). 
It suffices to verify the abstract assumptions (|4.7|)-(|4.9|): 



(4.7): If Wh\T G Pi, tlicn Vw/i is constant and Ih\'^Wh\ = |Vti;/i|. If Wh\T € Qi 
instead, then its gradient is linear and, consequently, jVui^l is convex. Then, 
we only need to realize that if the function ip is convex, then Xh^\T > ^\t- 
Indeed, using that Yl^=i ^j,T ^ 

A/V / aAt 

Ih^filriz) = V(p(zj,t)Aj,t(^) > <y3 y2zj,TXj^T{z) I = (p{z) 
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(4.8): (2.1), together with the L°°-stabihty of Uh, gives 

WUhWe - W\\l2 < \\TlhWe - w\\Ll\\IlhWe - w||loo < c\\w\\L'^\\IlhWe ~ w\\li. 

To bound the norm in the above inequaUty we add and subtract to 
obtain 

\\IlhW^ - w\\li < WUhWe - Wellii + \\We ~ w\\ < C ( — + £ ) |Dw|(il), 



where we have used (2.8), for p = 1, in conjunction with Proposition 2.1 



(4.9): We begin the proof of the energy diminishing property by 

■i'hinhw,) - *H = / ih\\/nhw,\ ~ \Dw\{n) < 

Jn 



I (i,,|vn^u.,| - |vn;,«;,|)') + ( I {\viihw,\ - \vw,\) 

|Vwe| - |Du;|(17)^ =1 + 11 + III. 



Applying Proposition 2.1 yields /// < ce|Dw|(il). We next invoke the triangle 
inequality along with interpolation estimate (2.8) for to arrive at 

f h 

II < \\7{nhW,-w,)\<ch\\D^w,\\Li <c-\Dw\{n). 
Jn e 

For the first term /, we use (2.8) to obtain 

/ = J^Ih\^UhW,\ - \^UhW,\ <chJ2 l|V(|Vn,,u;,|)||^i(j,) . 

For a smooth function, 
V(|V/|) = 



Ten 



^^D^f, a.e. {xeT: V/(x) ^ 0} , 



and V(|V/|) = a.e. {x e T : V/ 0}. This allows us to conclude 
that |V(|V/|)| < \D^f \ a.e. in T. This, together with the bound and 
Proposition |2.1[ shows that 



Jn e 



In 

The estimates above allow us to see that 

,2 \ 1/2 



£i(/i) = c he 



e2{h) = c 



Setting e = /i^/^ we obtain the result. □ 

Remark 4.7 (Energy diminishing interpolation). The proof of Corollary 
actually shows that if we were able to construct a TV diminishing interpolant, i.e., 
such that 

IvChM < \J^w\{n), 

then £2{h) = chfT^ and, setting e = h^^^ we would improve the rate of convergence to 
0{h^^^). Under the assumption that fl = (0, 1)'* and that the mesh is Cartesian, 
presents such a construction. 
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5. An Iterative Scheme for the Fully Discrete Scheme. A practical al- 
gorithm for the solution of the discrete variational inequalities ( |4.12 ) would require 
identifying the elements of the subdifferential d'i!h{u). In the continuous setting, re- 
lying on the results of {S^, such identification is presented in Let us describe this 
identification without going into technical details. We introduce the space 

X = {q e : Vq e L^{n), q-n = 0} , 

and stress that w E d'^{u) if and only if [2] 



and 



3zeXn«8i(L°°(l})) : w = -Vz. 



Notice that we can replace the first equality above by = z Du [2]. This serves 
as motivation for the discrete energy ^'/j given in ( 4.13[ ). Indeed, 



AAt A/'r 

Ten j=i Ten 3=1 qe«''Hq|<i 

= sup {Wwh,qh)h= sup {{'Vwh,qh)h-'3<3iiyiH)i^h)} , 

where 3s denotes the indicator function of S, the discrete space X/j is defined as 

X,, = {q,, e L°°(r!) : qh\T eQfyr eTh} , llq^llx, = max max ||q,,(z,.T)ll , 

Ten j=i^j\fj, 

and the discrete inner product (•, •)^ is defined by the quadrature rule 

Ten j=i 

The latter induces the norm ||q/i||^j = (qh,q/i)/,,, which clearly implies 

||q/i||h < c||q,i||L2, Vqhe^h- (5.1) 
In this setting it is not difiicult to see that the fully discrete subgradient fiow 



( |4.12[ ) with energy ( |4.13[ ) is equivalent to finding u^* C Yh and z^* C *Bi(X/i) that 
solve: 



(5.2) 



and 



(q„ - z;;+l, V<+1)^ < 0, Vq;, G 'Bi(X;,) 



(5.3) 



5.1. Exact Solver. Problem (^5.2h-(5.3) is not a practical numerical scheme 



since it involves the solution of the (local) variational inequality (5.3). To overcome 



this difiiculty we exploit that (|5.2|)-(|5.3|) are the optimality conditions of the functional 

fn- 



Wh t— 2St ll^^h ~ '^^llis + which motivates the following algorithm for the 
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solution of (5.2)-(5.3) (cf. [IllIS]): Let a > and r > be parameters to be chosen. 
Given g,, = < € V^, set = gt. A° = and find {vl,\l} C V,, x <Bi(X,0 by 



--5\ 



<0 Vq,. eSilXO, 



where v^^~^^ is the extrapolation defined in (2.3), and 



(5.4) 



(5.5) 



Finally, set = w,^. 

The last assertion needs to be justified, this is given in the following. 

Theorem 5.1 (Convergence of the inner iterations). Assume the inverse inequal- 
ity ||Vw;i||i2 < Ci/i~"'^|j?x;;i|ji2 holds, and the parameters t and a are chosen so that 
T < Cidh. Then the inner iterative loop (5.4) -(5.5 I converges, in the sense that for 
every L > 1, 



E 



L2 



< c. 



Proof. We follow the arguments of [H Proposition 3.1]. To alleviate the notation. 



set ei = ul+^ - vl and E'^ 



k+l 



A^. Subtract (5.5) from (5.2 1 to obtain 



Set q^ = "i^^^ in (5.4), q/j = Xj^^ in (5.3 1 and add them, to obtain 



-(5E' 



i+i 



V u' 



k+l 



V 



-pl + l 



< 0. 



(5.6) 



(5.7) 



Set Wh = '^Tc]^^ in (5.6) to obtain 



XII ' + 1||2 



o' + l||2 



i2+2r(E;+\Ve'+^) 



Multiply (|5.7|) by 2t and add it to (|5.8|) and add over ^ = 0, L — 1 to obtain 

2r , 



lie 



L||2 



-CrllE 



L\\2 
h Wh 



L-1 

E 



p/+l||2 



|H+^||i2+a||^E,-||K + — lie. 



;+iii2 



L-l 



(5.8) 



< e 



0||2 



-a||E0||^, + 2r5:(<52Ver,E'+\, 



Now we proceed as in [71 Proposition 3.1] with the last term: we use (2.5) to sum by 
parts and next employ inequality (5.1), the inverse inequality and repeated applica- 
tions of Cauchy-Schwarz to obtain the desired result. □ 

Remark 5.2 (Convergence of the Lagrange multiplier). As already mentioned in 
^ Remark 3.2(i)], convergence A', — > zj^^^^ . as I — oo, cannot he expected in general. 



Discrete TV Flows 



13 



due to the non-uniqueness o/z^^^. However the increments S\[^ converge to zero. 

Remark 5.3 (Solution of the discrete variational inequality). It is not difficult 
to see that the solution to the discrete variational inequality ( |5.4[ ) is explicit 



c{l,|AUfV<''+i|(z,,T)} 



5.2. Inexact Solver. It is not feasible to iterate in ( 5.4 )-( 5.5 1 until convergence. 



For this reason we included in the analysis presented in \ 4.2 a perturbation term since 



if we were to stop after L — 1 iterations, (5.5) would become 



1 

At 



1 



and (5.4) would read 



<0, 



with 

~k+i 



S^Wvf^. Doing the formal replacements vf^ 4— 



on the left hand side of these identities, setting Wh 



~k+l 



them we obtain 



■fe+i 

h 

Wh and adding 



and A^j 



8^ 
At 



*.(2|:+^)-*.(i5,)<(-^ 



7fe+l 



(5.9) 



ith 



*/i(w/i) 



Ten 



The non-homogeneous term on the right hand side can be understood as a pertur- 
bation p^*. In addition, the modified discrete energies can be controlled owing to 



Theorem 5.1 We make these ideas rigorous in the following. 

Theorem 5.4 (Convergence of fully discrete scheme with inexact solutions). Let 
he star-shaped with respect to a point and assume that Uq G BV{^) H L°°{^). If 
{?I^*,z^*} C V/i X are approximations to the solution of (5.2)-( |573| , computed 



with algorithm (5.4) -(5.5) in such a way that, for every time step k, the inequalities 



\T-Hvk\\L2<cAt^'\ 



< cAt, 



(5.10) 



are satisfied, then the following error estimate holds 



<||uo-<|U. +c(Ati/2 + /iV6). 



The spatial rate of convergence becomes 0{h^^^) as soon as Q — (0, 1)'' and the mesh 
Th is Cartesian. 

Proof. Since Uq e BV{n) n L°°{n), the solution, u'^*, to converges, with 

order 0{At^^^), to the exact solution. This also guarantees, via Corollary 
convergence of m^*, solution of (4.12), to u'^* with order 0{h^^^). 



4.6 



the 
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To study the convergence of u^* 



Wh 



in (4.121 and Wh 



,fe+i 



to u^* we proceed as in Theorem 4.5 Set 



we obtain 



in (5.9) and add the result. On denoting 



1 



Since 



< c r 



{Ih\^Wh + rj^\-Ih\^Wh\) < 



hWL^ 



Ten 



we notice that condition (5.10) suffices to guarantee that ||e 



< cM^I'^. Con- 



clude with a trivial application of the triangle inequality. 

The improvement on the spatial rate of convergence is due to Remark |4.7[ □ 



6. Numerical Experiments. To illustrate the theory developed in the preced- 
ing sections, here we present a series of numerical experiments. We implemented 



scheme (5.4)-(5.5) where the stopping criterion for the inner iterations is given by 



(5.10). The implementation was done with the help of the deal. II library (see 



IJISj). Unless noted otherwise, we set At = ^3^, r = At and a = 0.1 
6.1. The Characteristic of a Convex Set. Let (e C I 



where E is 



convex, connected and of finite perimeter. Define \e 



\E\ 



, where P{E) stands for 



the perimeter of E. Assume, in addition, that dE S C^'^ and that the curvature of 
E, X satisfies ||>f||L~(a£;) < ^e- Then, according to [211], if uq = Xe the solution to 



(3.2) with homogeneous Dirichlet boundary conditions is u{x,t) = (1 — A^t) XEix). 



Fig. 6.1. L 

solution u{x,t) = (1 



{L^)-errors for the evolution of the characteristi c yp of a circle B with exact 
\Bt)+XB{x) and\B = P{B)/\B\, fUgF (see^l ). 



Figure 6.1 



shows the L°°(L^)-error between the discrete and numerical solutions 
in the case when 17 — (—3, 3)^ and E — B{0, 1). As we can see the behavior of the 
error is actually better (i.e., 0{h^^'^)) than what our theory predicts. 



6.2. The Characteristic of a Disconnected Set. With the setting of §6.1[ 
references [H [8j also show that if Uq — J2'i^iXB{xi,ri)j where the centers are at the 
vertices of an equilateral triangle of unit side-length and the radii satisfy < 0.2, 
then u{x,t) = I]Li(l ~ ^B,t)+XB{x,,ri)- 




Fig. 6.2. Evolution of the characteristic set of three balls Bi = B{xi,ri) with centers Xi 
at the vertices of an equilateral triangle of unit size and radii ri < 0.2. The exact solution is 
u{x,t) = ~ ^Bit)'^XB{xi.ri) ■ T^^i-^ mesh is uniform with size h = and the time-step 

At = \/2/i/10. The discrete solution is shown every three time steps until extinction and preserves 
the structure of u without numerical diffusion and/or oscillations. 



Figure 6.2 shows the evolution in this case when h = 0.03. Notice that our method 
provides a good approximation with relatively few nodes and it does not introduce 
spurious oscillations. In addition, we obtain fairly good agreement with the extinction 
time, given by = = | = ^^j. In contrast to methods that involve regularization 
pn [55] , we observe consistency with the PDE in that the support of the fully discrete 
solution remains consta nt in time. Finally, the L°°(L^)-norm of the error with respect 
to h is shown in Figure 6.3 Again, we see that the error is 0{h^/'^). 



6.3. The Characteristic of an Annulus. Again in the setting of |6.1[ the 
evolution of Uq — Mxe, where E = B{0, R) \ B{0, r) is an annulus and M, i?, r S M, 
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Fig. 6.3. L°°{L^) -errors for the evolution of the characteristic of three circles (see The 
observed convergence rate Oih}^'^) is better that predicted by the theory. 



< R, is governed by [21 18] 



u(a;, t) 



'sgn(M)(|Af| - XEt) + XE + XB{0.r)tXB(0.r), t<Ti, 
^sgn(rn) (|m| - \B(Q,r){t - ^i))^ XB(o,fl), t > Ti, 



with Ti = 



\M\ 



and 



^B(0,r)Tl 



We set M = 4, R = 1/2, r = 1/4. Figure 6.4 shows the evolution of the numerical 
solution for h = 2^^ and notice that, in this case, Ti = 1/4, to = 2 and the extinction 
time is Text = 3/4. Numerically we obtained Ti^h ~ 0.2558, rrih ~ 2.045, Text,h ~ 
0.796; which are in good agreement with the exact values. In addition, we see that 
there is no diffusion, as opposed to what methods base d on regularization obtain. The 

Again we obtain 0{h^/'^). 



L°°{L )-eTroT with respect to h is presented in Figure 



6.5 



6.4. Convergence of the Inner Loop. The result of Theorem |5 . 1 1 guarantees 
convergence of the inner iteration at every step. However, there is no assessment of 
the speed of convergence. We numerically investigate the effect of the choice of initial 
condition for the inner loop. Set ~ (—2,2) and Ug = We choose 



4 = 0, 



x + 2, xe(-2,-l), 
-X, xG(-1,1), 



X 



2, xG(l,2). 



It is possible to show that only Zg e (95'(uo). For h — 2 ^, At = t = h/10 and 



cr = 1, Figure 6.6 plots the number of iterations as a function of the step for zj and 



Zg. The number of initial iteration is much higher (more than 5000) for Zq than for Zg. 
Moreover, notice the spike on both graphs at 79 steps. This is due to the fact that, 
at this step, extinction occurs. From this we conclude that the number of iterations 
heavily depends on the initial choice of zg or, more generally, on A^. Strategies for 
choosing it need further investigation. 

6.5. Comparison with Regularized Flow. To conclude our discussion, it is 
imperative to make a comparison between our method and those that involve regular- 
ization and show that the advantages in our approach are numerous. Before embarking 
in such an endeavor, let us recall some properties of regularized flows. The analysis 
developed in |22j does not provide a clear understanding of the relation between the 
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Fig. 6.4. Evolution of the characteristic of a ring E = B{0,R) \ B{0,r) with R = 1/2 and 
r = 1/4. The mesh is uniform with size h = and the time-step At = \/2h/W. The discrete 
solution, shown after 0, 18, 26, 50, 52, 72, 90, 108 and 170 time steps, reproduces well the form of 
the exact solution without numerical diffusion. The sealed characteristic Xb{o,r) decreases whereas 
Xs(o,r) increases until their heights merge and Xs(o,-R) continues to decrease. 

regularization e and the space discretization h. By means of numerical experiments 
the authors conclude that e = 0{h^) is the optimal scaling law, see [2^ Figures 4-5]. 

Let us, with the help of the theory developed i n §4.2[ try to bring some light into 
this matter. If we denote by the solution to (1.1) with regularization 5'e(w) = 
/j^ _^ |Vw|2 given by ([O]), [351, Theorem 2] shows that 



(6.1) 



is the solution to a fully discrete approximation of the regularized flow 



with "^ch = under the assumption that e i°°(0, T\ H^{Q)) we have 

||u, - utlh-iL^) < c i^t'^' + h'/h-'^/') . (6.2) 
To see this it suffices to realize that (4.7 1 is trivially satisfied and, setting Ch — H^, 



we conclude that (4.8) amounts to 

W^hW - w\\l2 < ch'^\\w\\H2 

and, since 



=^ £i(/i) = c/i^||t(;||^2, 
Jn Jn 
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Fig. 6.5. L°° (L^)- errors for the evolution of the characteristic of a ring (see 
experimental convergence rate 0{h^^^) is better than predicted. 



\6.3V. The 




Fig. 6.6. Number of iterations. Left: the initial dual variable is not in the subdifferential of the 
initial data, which requires more than 5000 initial iterations for convergence. Right: the ini tial d ual 
variable is in the subdifferential, which entails fewer initial iterations for convergence (see ^6.4^. 



we have, for (4.9 ), 



£2{h) < ch\\w\\w^. 

Finally, [5TJ Theorem 1.2] shows that ||ue||i2(jLf2) < e^" for some a £ Nq. Combining 
(6.1) and (6.2) with Theorem 4.5 we conclude 



ufX\\L--{L^) < C 



(A<1/2 + el/2 ^ ^ ^1/2 



(6.3) 



which yields the optimal scaling e = ©(/it+o). 

The last ingredient we need to make the comparison is to recall that (cf. pSl 
Remark 5.3]) in the one dimensional case, i.e., d = 1, if the initial data is monotone, 
it is itself a minimizer of the total variation energy, and so the flow fixes it. In other 
words, if Uq is monotone, then u{t) = Uq for all t > 0. 

Consider, in f2 = (0, 1), the initial data 



Uo(a;) 



i<x<l. 



According to the discussion presented above, the solution to (3.2) is u{t) = Uq. Fig- 
ure 6.7 shows the solution, at T = 5, obtained with our method and the regularized 
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0.2 0.4 0.6 o.a 1 



Fig. 6.7. Solution, at T = 5, of the TV flow with monotone initial data. Our method (red 
with +), regularized flow with e = h (blue with X ) and with e = h^ (magenta with *). h = 2~^, 
At = 2"-'''. The unregularized solution coincides with the exact solution, whereas the regularized 
flow misses it, although it belongs to the discrete space. 



flow with € = 0{h ) and e = 0{h). The first choice is the one advocated in [55]; while 



Proposition 2.1 In such a case, (6.3) gives the foflowing error estimate 



the second is the optimal according to (6.3) provided a = 1, which is consistent with 



provided e k, h ~ At. The mesh size is h = and the time-step At = 2~^°. 
Notice that the requirement At = 0{h^), which is needed for L^-convergence in the 
regularized flow (cf. [551 Theorem 4] and [5TJ Theorem 1.7]), is satisfied. 

The advantages of our method are now evident. We do not impose any restric- 
tion on the time-step, as opposed to the At = 0{h?) that is necessary in regularized 
methods to guarantee convergence in L^(J7). Even if one is willing to settle for conver- 
gence in LP(f2), with p < d/ld— 1), the methods with regularization require that the 
solution of the regularized flow is in L°° (0, T; Wl (n)) n (0, T; Hl^^{n)). If an error 
estimate is desired, one must impose that Uq G C^(ri), and even in that case it is not 
clear what is the relation between /i. At and e. In addition to these approximation 
issues, the regularized flow requires the solution, at each time step, of a nonlinear 
system and no convergence analysis is provided. In contrast, we have developed and 
analyzed an inexact iterative scheme for the solution of our problems at each time 
step, and we have showed its global convergence. To finalize, the result presented in 
Figure |6.7| shows that the regularized flow misses certain fundamental features of the 
problem, even in simple cases. 

7. Total Variation Minimization. We conclude with yet another application 



of our result on approximation of functions of bounded variation (Proposition 2.1 ): we 
improve on the existing results about total variation minimization. Let g e i°°(f2) 
and a > 0. Consider 'B.{w) — ^{w) + ^\\w — (?||^2. Thanks to the fact that this 
functional is strictly convex there exists a unique ^ G BV{il.) L'^iU) such that 

S(^) = inf {S(w) : w e L^{VL)] < oo. (7.1) 

Here we are interested in the approximation of ^ by elements of V^- Since V/i is finite 
dimensional, there is a unique G V/i such that 

S(60 = inf {^{wh) : in e V4 < oo. (7.2) 



20 



S. Bartels, R.H. Nochetto, A.J. Salgado 



The main approximation properties of are detailed in the following. 

Theorem 7.1 (Convergence of discrete minimizers). Assume that fl is star 
shaped with respect to a point. Let ^ € BV(fl) D L°°(D,) and G Wh be defined as in 
(7.1) and (7.2 1, respectively. Then 

U-^h\\L^<ch^/^. 



Proof. We adapt the ideas presented in [71 Theorem 3.1]. Let e > and 
£,e G C°°(f2) be an approximation of ^ that satisfies all the properties stated in Propo- 
sition |2.1| Owing to the strict convexity of S and the fact that ^/j is a discrete 
minimizer, we have 

f lie - aili. < - m) < s(n,c.) - m 

= (II vn.CelUi - mm) + f (lin.Ce - 5lli. - U gWh) = + ^2, 

where Hh is the Clement interpolation operator |15 . Let us look at each one of the 
terms in this last inequality: 

Ai'. We add and subtract the W^;^^-seminorm of and use its approximation properties 
with respect to total variation along with the approximation properties of 11^ 
described in (2.8 1: 

Ai < \\s/{n^^, - 6)iili + iive.iiLi - mm) < c + mm)- 

A2- We, again, use the approximation properties of and the fact that ^ and g are 
essentially bounded, say by a constant c > 0, 

^2 < f (III.6 - e.iui + lie. - eiuo <c(^^+e^ mm)- 

Setting e = h^ we obtain the result. □ 

Remark 7.2 (Convergence of total variation minimization). Theorem \ 7. 1\ is, 
in a sense, an improvement over the original result of 171, at least for star .shaped 
domains, and under the boundedness assumptions on g and If S, & B^{L'^{fl)) for 
some s e (0, 1], and relying on the results of \S7^ , ^ Theorem 3.1] proves the estimate 

U-^hWv^ < ch^, 

so that the best possible rate of convergence is 0{h^^'^), which is what we obtain, but 
with lower regularity. To understand this regularity assumption it suffices to recall 
that BV{Vl) n L°°(rj) B^{L^{n)) for s < ^ (see JM Lemma 38.1] for a proof 
and, in some sense, the converse inclusion). The key step that allowed us to reduce 
the smoothness assumption is Proposition \2. 1\ In addition, the proof of Theorem \7.1\ 
shows that if we had a TV-diminishing interpolant, we would obtain 

U'ih\\L-<ch^'^ 

which is an optimal error estimate for ^ G i?^(L^(f2)). Such a construction is pre- 
sented in l3£ 
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